#ifndef ONEDIMDAMBREAK_HPP
#define ONEDIMDAMBREAK_HPP

#include "edgeTraversal.hpp"
#include "sweSolver.hpp"
#include "solverBase.hpp"
#include "meshTypes.hpp"
#include "meshFunctions.hpp"
#include "magic_enum.hpp"
#include "vtkFunctions_export.hpp"

#include "solverSettings.hpp"


#include <filesystem>
namespace fs = std::filesystem;

using namespace hms;

struct oneDimDamBreak{

	hms::Time time;
	hms::Traversal traversal;
	hms::EdgeTraversal edgeTraversal;
	hms::CellTraversal cellTraversal;
	Index blocksize;
	std::string fluxSchemeStr;
	bool
		useMuscl,
		useExport,
		useAtomicWrite,
		dryCase;
	scalar maxCourant;
	fs::path exportPath;

	oneDimDamBreak( hms::Time&& t, hms::SolverSettings&& s, bool dryCaseArg ) :
		time { std::move(t) },
		traversal { s.traversal },
		edgeTraversal { s.edgeTraversal },
		cellTraversal { s.cellTraversal },
		blocksize { s.blocksize },
		fluxSchemeStr { std::move(s.fluxScheme) },
		// useFriction { s.useFriction },
		// useSourceFields { s.useSourceFields },
		useMuscl { s.useMuscl },
		useExport { s.useExport },
		useAtomicWrite { s.useAtomicWrite },
		dryCase { dryCaseArg },
		maxCourant { s.maxCourant },
		exportPath { std::move(s.exportPath) }
	{}

	template<typename MeshType>
	void operator()(MeshType& mesh) {
	
		/**************************
		***** case properties *****
		***************************/
		bool
			useFriction     {false},
			useSourceFields {false};
		scalar 
			dLeft           {4.0},
			dRight          {(dryCase) ? 0.0 : 1.0};
		
		LimiterFunction lf  {LimiterFunction::minMod};

		static_assert( std::is_same_v<MeshType, UniMesh> || std:: is_same_v<MeshType, RectMesh>, 
		"1D Dam Break Verification case can only run on MeshTypes: RectMesh, UniMesh. The entered type does not match. Please choose on of the two options.");

		std::cout<< "Verification test case: 1D dam break " <<((dryCase) ? "(dry bed)" : "(wet bed)") << std::endl; 

		

		/**************************
		*********** mesh **********
		***************************/
		Index
			nx { mesh.nCellsX() },
			ny { mesh.nCellsY() },
			ng { mesh.nGhostLayers() };

		mesh.blocksize(blocksize);

		static_assert( std::is_same_v<MeshType, UniMesh> || std:: is_same_v<MeshType, RectMesh>, 
		"Rainfall Runoff Verification case can only run on MeshTypes: RectMesh, UniMesh. The entered type does not match. Please choose on of the two options.");
		
		Index
			cellsPerRow {nx + 2*ng},
			cellsPerCol {ny + 2*ng};

		/**************************
		*********** Time **********
		***************************/
		// Time time;

		// time.startAt(0.0);
		// time.endAt(simulationTime);
		// time.maxTimeStep(maxTimeStep);
		// time.writeAtTimeIntervals(useWriteTimeInterval);
		// if ( useWriteTimeInterval ) time.writeTimeInterval(writeInterval);
		// else time.writeStepInterval(static_cast<hms::Index>(writeInterval));

		/**************************
		**** solver properties ****
		***************************/
		SWESolver<MeshType> solver {mesh, time};

		using FluxScheme = typename SolverTraits<SWESolver<MeshType>>::FluxScheme;
		FluxScheme fluxScheme = FluxScheme::HLLC;

		if(!fluxSchemeStr.empty()) {
			auto fluxSchemeOpt { magic_enum::enum_cast<FluxScheme>(fluxSchemeStr) };
			if ( fluxSchemeOpt.has_value() )
			fluxScheme = fluxSchemeOpt.value();
		}
		else {
			std::cerr << "WARNING: Flux Scheme \"" << fluxSchemeStr
				<< "\" unknown! Defaulting to HLLC\n";
			fluxScheme = FluxScheme::HLLC;
		}

		solver.useAtomicWrite(useAtomicWrite);
		solver.time(time);
		solver.fluxScheme(fluxScheme);
		solver.traversal(traversal);
		solver.edgeTraversal(edgeTraversal);
		solver.cellTraversal(cellTraversal);
		solver.useFriction(useFriction);
		solver.useSourceFields(useSourceFields);
		solver.useMuscl(useMuscl);
		solver.limiterFunction(lf);
		solver.maxCourant(maxCourant);

		/**************************
		******** fields ***********
		***************************/
		solver.v().values() = 0.0;
		solver.z().values() = 0.0;

		solver.d().values().reshaped(cellsPerRow, cellsPerCol).
			transpose().rightCols(cellsPerRow/2) = dRight;
		solver.d().values().reshaped(cellsPerRow, cellsPerCol).
			transpose().leftCols(cellsPerRow/2)  = dLeft;

		/**************************
		**** boundary patches *****
		***************************/
		const BoundaryPatch& left   { mesh.boundaryPatch("left") };
		const BoundaryPatch& right  { mesh.boundaryPatch("right") };
		const BoundaryPatch& bottom { mesh.boundaryPatch("bottom") };
		const BoundaryPatch& top    { mesh.boundaryPatch("top") };

		/**************************
		*** boundary conditions ***
		***************************/
		// bottom elevation
		
		solver.z().setBoundaryCondition( left,   StaticValue{}, Array1Ns<1>{0}, true );
		solver.z().setBoundaryCondition( right,  StaticValue{}, Array1Ns<1>{0}, true );
		solver.z().setBoundaryCondition( top,    StaticValue{}, Array1Ns<1>{0}, true );
		solver.z().setBoundaryCondition( bottom, StaticValue{}, Array1Ns<1>{0}, true );
		// water depth
		solver.d().setBoundaryCondition( left,   ZeroGradient{} );
		solver.d().setBoundaryCondition( right,  ZeroGradient{} );
		solver.d().setBoundaryCondition( top,    ZeroGradient{} );
		solver.d().setBoundaryCondition( bottom, ZeroGradient{} );
		// flow velocity
		solver.v().setBoundaryCondition( left,   Slip_velocity{} );
		solver.v().setBoundaryCondition( right,  Slip_velocity{} );
		solver.v().setBoundaryCondition( top,    Slip_velocity{} );
		solver.v().setBoundaryCondition( bottom, Slip_velocity{} );
		

		/**************************
		***** initialisation ******
		***************************/
		
		solver.initialise(); 
		mesh.printInfo();
		std::cout << "\nFriction is enabled: " << std::boolalpha << useFriction
			<< "\nSource fields are enabled: " << useSourceFields
			<< "\nUsing second order reconstruction: " << useMuscl
			<< "\n";

		if(traversal == hms::Traversal::edges) {
			std::cout << "Edge traversal: "
				<< magic_enum::enum_name( solver.edgeTraversal() ) << std::endl;
		} else {
			std::cout << "Cell traversal: "
				<< magic_enum::enum_name( solver.cellTraversal() ) << std::endl;
		}
		std::cout << "flux scheme: " << magic_enum::enum_name( solver.fluxScheme() ) <<std::endl;

		Timer t;
		std::cout << std::setprecision(16) << std::fixed;
		std::cout << "Time: " << time.current() << "s";

		std::unique_ptr<vtk::xml::XMLDocument> bindingFile {(useExport) ? vtk::createBindingDoc() : nullptr };

		while ( !time.atEnd() ){
			// vtk export
			if( useExport && time.atWriteTime()) {
				vtk::toVtk(exportPath,mesh,solver,time.nWritten());
				vtk::updateBinding(bindingFile.get(), mesh, time.current(), time.nWritten(), exportPath);
			}

			// evolution in time.
			solver.advance( time );
		}

		// last vtk export
		if( useExport ){
			vtk::toVtk(exportPath,mesh,solver,time.nWritten());
			vtk::updateBinding(bindingFile.get(), mesh, time.current(), time.nWritten(), exportPath);
		} 

		double elapsed { t.elapsed() };
		double execTimeSum {0};
		for ( auto execTime : solver.execTimings() )
			execTimeSum += execTime;

		std::cout
			<< "\nComputation time: " << elapsed << "s"
			<< "\nSum of execution times logged by solver: " << execTimeSum << "s"
			<< "\n\t" << 100 * solver.execTimings()[0] / execTimeSum
				<< "% for time step calculation"
				<< " (" << solver.execTimings()[0] << "s)"
			<< "\n\t" << 100 * solver.execTimings()[1] / execTimeSum
				<< "% for fluxes (including slope)"
				<< " (" << solver.execTimings()[1] << "s)"
			<< "\n\t" << 100 * solver.execTimings()[2] / execTimeSum
				<< "% for sources " << ( useFriction ?
					std::string{ "(including friction)" } :
					std::string{ "(friction disabled)" } )
				<< " (" << solver.execTimings()[2] << "s)"	
			<< "\n\t" << 100 * solver.execTimings()[3] / execTimeSum
				<< "% for new state"
				<< " (" << solver.execTimings()[3] << "s)"
			<< "\n\t" << 100 * solver.execTimings()[4] / execTimeSum
				<< "% for updating boundaries"
				<< " (" << solver.execTimings()[4] << "s)"
			<< "\n";
	}
};

#endif